      subroutine ornl_dcor2d(    f,incf,ldf,ifx0,n_fx,ify0,n_fy,
     $			    g,incg,ldg,igx0,n_gx,igy0,n_gy,
     $			    h,inch,ldh,ihx0,n_hx,ihy0,n_hy)
      double precision f(ldf,n_fy)
      double precision g(ldg,n_gy)
      double precision h(ldh,n_hy)
      integer incf,ldf,ifx0,n_fx,ify0,n_fy
      integer incg,ldg,igx0,n_gx,igy0,n_gy
      integer inch,ldh,ihx0,n_hx,ihy0,n_hy
c-----------------------------------------------------------------------------
c
c   dcor2d  performs a 2D correlation in the Space domain :
c	h(i,j) = Sum[ f(k,l) * g(i+k,j+l) ]
c
c-----------------------------------------------------------------------------
c
c   PARAMETERS:
c
c	f:	Vector containing sequence "f"
c	ldf:	Increment between two successive columns of "f"
c	nfx1:	Index of the first element of each column of "f"
c	nfx2:	Index of the  last element of each column of "f"
c	nfy1:	Index of the first column of "f"
c	nfy2:	Index of the  last column of "f"
c	
c	g:	Vector containing sequence "g"
c	ldg:	Increment between two successive columns of "g"
c	ngx1:	Index of the first element of each column of "g"
c	ngx2:	Index of the  last element of each column of "g"
c	ngy1:	Index of the first column of "g"
c	ngy2:	Index of the  last column of "g"
c	
c	h:	Vector containing sequence "h"
c	ldh:	Increment between two successive columns of "h"
c	nhx1:	Index of the first element of each column of "h"
c	nhx2:	Index of the  last element of each column of "h"
c	nhy1:	Index of the first column of "h"
c	nhy2:	Index of the  last column of "h"
c	
c-----------------------------------------------------------------------------
c
c   PLEASE NOTE:
c	Please note that the array pointers must all point to the first 
c	element of the array "(nfx1,nfy1)", "(ngx1,ngy1)" and "(nhx1,nhy1)"
c	 If "f" for example is defined as
c	    dimension f(-25:45,10:21)
c	Then "dcor2d" must be called with the following parameters
c	    call dcor2d( f(-25,10),(45-(-25)+1),-25,45,10,21 ... )
c
c-----------------------------------------------------------------------------
c
c   This Fortran Subroutine written by
c	Jean-Pierre Panziera
c	Silicon Graphics Inc
c	September 9, 1991
c
c-----------------------------------------------------------------------------

      double precision tmp, zero, one
      parameter (zero = 0.0, one = 1.0)
      integer i,j,k,j0,j1,k1,k2,ify1,igy1,ihy1,ix,jh,jf,jg,ig,igx1
c-----------------------------------------------------------------------------
c
c   Compute Range
c
c-----------------------------------------------------------------------------
	ify1 = ify0 + n_fy - 1
	igy1 = igy0 + n_gy - 1
	ihy1 = ihy0 + n_hy - 1

	igx1 = igx0 + n_gx - 1
	ig   = 1 + (igx1-igx0) * incg

	j0 = max( ihy0, ify0-igy1)
	j1 = min( ihy1, ify1-igy0)
c-----------------------------------------------------------------------------
c
c   Zero the Resulting sequence
c
c-----------------------------------------------------------------------------
	do j = ihy0, ihy1
	    ix = 1
	    jh = j-ihy0 + 1
	    do i = 1, n_hx
	        h(ix,jh) = zero
		ix = ix + inch
	    end do
	end do
c-----------------------------------------------------------------------------
c
c   Then Compute the correlation
c
c-----------------------------------------------------------------------------
	do j = j0, j1
	    k1 = max( -igy1, j-ify1)
	    k2 = min( j-ify0, -igy0)
	    jh = j -ihy0+1
	    do k = k1, k2
		jf = j-k -ify0+1
		jg = -k  -igy0+1
		call dfir1d(	f(1,jf), incf, ifx0, n_fx,
     $				g(ig,jg), -incg, -igx1, n_gx,
     $				h(1,jh), inch, ihx0, n_hx,
     $				one, one)
	    end do
	end do
c-----------------------------------------------------------------------------
      return
      end
